function QBFmat = BeamformingSubarray(qbit,L,TRX,fc,BW_des,Ptx,showfig,qSubarray)



x1 = TRX(1);y1 = TRX(2);z1 = TRX(3);
x2 = TRX(4);y2 = TRX(5);z2 = TRX(6);
z12 = z1^2;z22 = z2^2;
% 
% fc = 30e9;
c = 3e8;
% BW_des = 1.5e9;
% Ptx = 0;%dBm


Ptx_lin = 10.^((Ptx-30)/10); % in watts
Lb = c/fc;
spac = Lb/2;
hp = 6.625e-34;
kb = 1.38e-23;
Tk = 290;


LIS_f = 0;
N = round(L/spac);

xcoord=linspace(-L/2,L/2,N);
ycoord=linspace(-L/2,L/2,N);
[xx,yy]=meshgrid(xcoord,ycoord);
coord_mat=xx+1i*yy;


a = (sqrt((xx-x1).^2+(yy-y1).^2+z1^2)+sqrt((xx-x2).^2+(yy-y2).^2+z2^2))/2;

Distmat = a*2;



% Hmat=(Lb/(2*pi))*exp(-1i*2*pi*fc*Distmat/c)./(Distmat1)./Distmat2;% channel response
% %     Hmat(Dist_arr>Radius)=0;
% Nactiv=size(a(:),1);%valid number of points
% BFmat=exp(-1i*angle(Hmat));% traditional angle design
if length(qSubarray) == 1
    fsub = linspace(fc-BW_des/2,fc+BW_des/2,2^qSubarray);
    fsub = repmat(fsub,N^2/2^qSubarray,1);
    fsub = reshape(fsub,N,N);
else 
    fsub = linspace(fc-BW_des/2,fc+BW_des/2,2^sum(qSubarray));
    fsub = reshape(fsub,2^qSubarray(1),2^qSubarray(2));
    fsub = kron(fsub,ones(N/2^qSubarray(1),N/2^qSubarray(2)));
    % fsub = reshape(fsub,N,N);
end
BFmat = exp(-2j*pi*fsub.*Distmat/c);

QBFmat=quantz(exp(-1i*angle(BFmat)),qbit);% beamforming matrix


























end